%%%%Parametric Bootstrap of Tobit Parameters
rng(round(pi*1000000))

S          = 10000000 ; %%%%This is the number of parametric bootstrap samples
betaregPBS = nan(S,length(betareg)) ; %%%%This will hold parametric bootstrap draws of the Tobit Parameter Vector
Zsamp      = randn(S,K_E+K_L+12) ; %%%This is a matrix of standard normal draws, where each row will correspond to a simulated betareg vector
Vchol      = chol(Vhat) ;
goodsim    = zeros(S,1) ;
for ss=1:S
    betasamp = ( betareg' + Zsamp(ss,:)*Vchol ) ;
    betaregPBS(ss,:) = betasamp ;
    test1 = [betasamp(K_E+K_L+1); 
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+2)); 
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+3)); 
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+4));
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+2)+betasamp(K_E+K_L+3));
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+2)+betasamp(K_E+K_L+4));
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+3)+betasamp(K_E+K_L+4));
     (betasamp(K_E+K_L+1)+betasamp(K_E+K_L+2)+betasamp(K_E+K_L+3)+betasamp(K_E+K_L+4))] ;
    testflag1 = min(test1)<=0 ;
    test2 = [betasamp(K_E+K_L+5); 
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+6)); 
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+7)); 
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+8));
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+6)+betasamp(K_E+K_L+7));
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+6)+betasamp(K_E+K_L+8));
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+7)+betasamp(K_E+K_L+8));
     (betasamp(K_E+K_L+5)+betasamp(K_E+K_L+6)+betasamp(K_E+K_L+7)+betasamp(K_E+K_L+8))] ;
    testflag2 = min(test2)<=0 ;
    test3 = [betasamp(K_E+K_L+9); 
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+10)); 
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+11)); 
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+12));
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+10)+betasamp(K_E+K_L+11));
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+10)+betasamp(K_E+K_L+12));
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+11)+betasamp(K_E+K_L+12));
     (betasamp(K_E+K_L+9)+betasamp(K_E+K_L+10)+betasamp(K_E+K_L+11)+betasamp(K_E+K_L+12))] ;
    testflag3 = min(test3)<-1 | max(test3)>1 ;

    if (testflag1+testflag2+testflag3)==0
        goodsim(ss) = 1 ;
    end

% % % % % % % % % %     ss
% % % % % % % % % %     if (testflag1+testflag2+testflag3)==0
% % % % % % % % % %         disp('This one is good.')
% % % % % % % % % %     else
% % % % % % % % % %         disp('This one fails the positive definite condition.......................................................')
% % % % % % % % % %         if testflag1==1
% % % % % % % % % %             disp('NEGATIVE log(THETA_E) VARIANCES GENERATED.')
% % % % % % % % % %         end
% % % % % % % % % %         if testflag2==1
% % % % % % % % % %             disp('NEGATIVE log(THETA_L) VARIANCES GENERATED.')
% % % % % % % % % %         end
% % % % % % % % % %         if testflag3==1
% % % % % % % % % %             disp('INFEASIBLE CORRELATIONS GENERATED.')
% % % % % % % % % %         end
% % % % % % % % % %     end
% % % % % % % % % %     pause
end

format shortg
clear Zsamp
sum(goodsim)
betaregPBS(goodsim==0,:) = [] ;
% % % % % % % % % % for ss=1:length(betareg) 
% % % % % % % % % %     figure; 
% % % % % % % % % %     histogram(betaregPBS(:,ss),50); 
% % % % % % % % % %     hold on; 
% % % % % % % % % %     plot([betareg(ss),betareg(ss)],[0,sum(goodsim)/100],'LineWidth',4); 
% % % % % % % % % %     xlabel(num2str(betareg(ss)-mean(betaregPBS(:,ss)))); 
% % % % % % % % % % end

[max(max(abs((betareg'-mean(betaregPBS))))) mean(mean(abs((betareg'-mean(betaregPBS)))))] 
[max(max(abs((Vhat-cov(betaregPBS))))) mean(mean(abs((Vhat-cov(betaregPBS)))))]  
[max(max(abs((diag(Vhat)-diag(cov(betaregPBS)))./diag(Vhat)))) max(max(abs((diag(Vhat)-diag(cov(betaregPBS)))))) mean(mean(abs((diag(Vhat)-diag(cov(betaregPBS))))))]  
% % % % % % % % % % [diag(Vhat) diag(cov(betaregPBS))] 

getrows         = randsample(length(betaregPBS(:,1)),NBS) ; 
tempbetaregsamp = betaregPBS(getrows,:) ; 
fitnessval      = max(abs(betareg' - mean(tempbetaregsamp))) + max(max(abs(Vhat - cov(tempbetaregsamp))))   % mean(abs(betareg' - mean(tempbetaregsamp))) + mean(mean(abs(Vhat - cov(tempbetaregsamp))))   % sum((betareg' - mean(tempbetaregsamp)).^2) + sum(sum( (Vhat - cov(tempbetaregsamp)).^2 )) 
for ss = 1:100000
    getrows2         = randsample(length(betaregPBS(:,1)),NBS) ; 
    tempbetaregsamp2 = betaregPBS(getrows2,:) ; 
    fitnessval2      = max(abs(betareg' - mean(tempbetaregsamp2))) + max(max(abs(Vhat - cov(tempbetaregsamp2))))  ; % mean(abs(betareg' - mean(tempbetaregsamp2))) + mean(mean(abs(Vhat - cov(tempbetaregsamp2))))  ; % sum((betareg' - mean(tempbetaregsamp2)).^2) + sum(sum( (Vhat - cov(tempbetaregsamp2)).^2 )) ; 
    if fitnessval2<fitnessval
        getrows = getrows2 ; 
        tempbetaregsamp = tempbetaregsamp2 ; 
        fitnessval      = fitnessval2 ; 
    end 
    if mod(ss,5000)==0
        [ss fitnessval]
    end
end 
fitnessval
[max(max(abs((betareg'-mean(tempbetaregsamp))))) mean(mean(abs((betareg'-mean(tempbetaregsamp)))))] 
[max(max(abs((Vhat-cov(tempbetaregsamp))))) mean(mean(abs((Vhat-cov(tempbetaregsamp)))))]  
[max(max(abs((diag(Vhat)-diag(cov(tempbetaregsamp)))./diag(Vhat)))) max(max(abs((diag(Vhat)-diag(cov(tempbetaregsamp)))))) mean(mean(abs((diag(Vhat)-diag(cov(tempbetaregsamp))))))]  


betaregPBS = tempbetaregsamp ;
clear tempbetaregsamp tempbetaregsamp2 getrows getrows2
goodsim    = length(betaregPBS(:,1)) ;
sum(goodsim)
for ss=1:length(betareg) 
    figure; 
    histogram(betaregPBS(:,ss),35); 
    hold on; 
    plot([betareg(ss),betareg(ss)],[0,sum(goodsim)/100],'LineWidth',4); 
    xlabel(num2str(betareg(ss)-mean(betaregPBS(:,ss)))); 
end

if Model==4.05
    save('.\STAGE4\betaregPBS_Model4point05.mat','betaregPBS')
elseif Model==4.085
    save('.\STAGE4\betaregPBS_Model4point085.mat','betaregPBS')
end